## BELL (2024) ESTIMATOR ON ISSP DATA

# Create Dummies
clean_ISSP_data_Bell <- fastDummies::dummy_cols(clean_ISSP_data, select_columns = c("wage_q", "country_s"))

# Group Controls
wage_controls <- colnames(clean_ISSP_data_Bell)[grep("^wage_q_", colnames(clean_ISSP_data_Bell))]
wage_controls <- wage_controls[!wage_controls %in% "wage_q_1"] # omitted category for regression

country_controls <- colnames(clean_ISSP_data_Bell)[grep("^country_s_", colnames(clean_ISSP_data_Bell))]
country_controls <- country_controls[!country_controls %in% "country_s_AT"] # omitted country for regression

controls <- paste(c(wage_controls, country_controls), collapse = " + ")
regression_equation_ISSP <- paste("education ~ 0 + high_meaning + high_flex + high_telecommute + ", controls)

# First Stage
first_stage <- lm(formula = regression_equation_ISSP, weights = weight, data = clean_ISSP_data_Bell)
rss_full <- sum(residuals(first_stage)^2)
n = nobs(first_stage)
k = qr(model.matrix(first_stage))$rank

# Save First Stage for Table
first_stage_issp <- first_stage

# Compute Coefficients + Anderson-Rubin bounds
pi <- coef(first_stage)["high_meaning"]
delta <- coef(first_stage)

cov_matrix <- vcov(first_stage)

c_pipi <- cov_matrix["high_meaning", "high_meaning"]
c_deltadelta <- diag(cov_matrix)
c_delpi <-  cov_matrix["high_meaning", ]

crit <- 1.96

a <- pi^2 - (crit^2 * c_pipi)
b <- 2 * crit^2 * c_delpi - 2*delta*pi
c <- ( delta^2 - (crit^2)*c_deltadelta)

lb <- (- (-b + sqrt(b^2 - 4 * a * c)) / (2 * a))
ub <- (- (-b - sqrt(b^2 - 4 * a * c)) / (2 * a))

beta <- -delta/pi
  
# Nested First Stage
regression_equation_nest <- paste("education ~ 0 + high_flex + high_telecommute +", controls)
first_stage_nest <- lm(formula = regression_equation_nest, weights = weight, data = clean_ISSP_data_Bell)
rss_nest <- sum(residuals(first_stage_nest)^2)
  
# Compute partial F 
partial_F = (rss_nest - rss_full) / (rss_full / (n - k))

# Linear Regressions (base and with productivity controls) + Bell Estimates
controls_country <- paste(c(country_controls), collapse = " + ")

base_equation <- paste("high_meaning ~ high_flex + high_telecommute +", controls_country)
base_ISSP <- lm(base_equation, 
                weights=weight,
                data = clean_ISSP_data_Bell)
  
prod_equation <- paste("high_meaning ~ high_flex + high_telecommute + education + ", controls_country)
prod_ISSP <- lm(prod_equation, weights=weight, data = clean_ISSP_data_Bell)

# Confidence intervals
conf_base_ISSP <- confint(base_ISSP, level = 0.95)
conf_prod_ISSP <- confint(prod_ISSP, level = 0.95)
  
# Create Table
base_flex_ISSP <- data.frame(conf_base_ISSP["high_flex",1], base_ISSP$coefficients["high_flex"], conf_base_ISSP["high_flex",2]) 
base_flex_ISSP[] <- lapply(base_flex_ISSP, function(x) {if (is.numeric(x)) {format(round(x, digits = 3), nsmall = 3)} else {x}})
  
base_tel_ISSP <- data.frame(conf_base_ISSP["high_telecommute",1], base_ISSP$coefficients["high_telecommute"], conf_base_ISSP["high_telecommute",2])
base_tel_ISSP[] <- lapply(base_tel_ISSP, function(x) {if (is.numeric(x)) {format(round(x, digits = 3), nsmall = 3)} else {x}})
  
prod_flex_ISSP <- data.frame(conf_prod_ISSP["high_flex",1], prod_ISSP$coefficients["high_flex"], conf_prod_ISSP["high_flex",2]) 
prod_flex_ISSP[] <- lapply(prod_flex_ISSP, function(x) {if (is.numeric(x)) {format(round(x, digits = 3), nsmall = 3)} else {x}})
  
prod_tel_ISSP <- data.frame(conf_prod_ISSP["high_telecommute",1], prod_ISSP$coefficients["high_telecommute"], conf_prod_ISSP["high_telecommute",2])
prod_tel_ISSP[] <- lapply(prod_tel_ISSP, function(x) {if (is.numeric(x)) {format(round(x, digits = 3), nsmall = 3)} else {x}})
  
Bell_flex_ISSP <- data.frame(lb["high_flex"], beta["high_flex"], ub["high_flex"])
Bell_flex_ISSP[] <- lapply(Bell_flex_ISSP, function(x) {if (is.numeric(x)) {format(round(x, digits = 3), nsmall = 3)} else {x}})
  
Bell_tel_ISSP <- data.frame(lb["high_telecommute"], beta["high_telecommute"], ub["high_telecommute"])
Bell_tel_ISSP[] <- lapply(Bell_tel_ISSP, function(x) {if (is.numeric(x)) {format(round(x, digits = 3), nsmall = 3)} else {x}})
  
  
## Open a connection to a file
file_conn <- file(file.path(graph_dir, "Proxy_ISSP.txt"), "w")
  
# Write the LaTeX code to the file
writeLines("\\begin{tabular}{llll}", file_conn)
writeLines("\\hline \\hline", file_conn)
writeLines("              & Base        & Productivity Controls & Bell Proxy  \\\\", file_conn)
writeLines("\\hline", file_conn)

# Write point estimates for TELECOMMUTE
writeLines(paste("Telecommuting", " & ",
                 base_tel_ISSP[1,2], " & ", 
                 prod_tel_ISSP[1,2], " & ", 
                 Bell_tel_ISSP[1,2], " \\\\"), file_conn)

# Write confidence intervals for TELECOMMUTE
writeLines(paste("   \\textit{Conf. Int.}           & ", 
                 "(", base_tel_ISSP[1,1], ",", base_tel_ISSP[1,3], ")", " & ",
                 "(", prod_tel_ISSP[1,1], ",", prod_tel_ISSP[1,3], ")", " & ",
                 "(", Bell_tel_ISSP[1,1], ",", Bell_tel_ISSP[1,3], ")", " \\\\[0.5em]"), file_conn)

# Write point estimates for FLEXIBILITY
writeLines(paste("High Flexibility", " & ",
                 base_flex_ISSP[1,2], " & ", 
                 prod_flex_ISSP[1,2], " & ", 
                 Bell_flex_ISSP[1,2], " \\\\"), file_conn)

# Write confidence intervals for FLEXIBILITY
writeLines(paste("   \\textit{Conf. Int.}           & ", 
                 "(", base_flex_ISSP[1,1], ",", base_flex_ISSP[1,3], ")", " & ",
                 "(", prod_flex_ISSP[1,1], ",", prod_flex_ISSP[1,3], ")", " & ",
                 "(", Bell_flex_ISSP[1,1], ",", Bell_flex_ISSP[1,3], ")", " \\\\[0.25em]"), file_conn)

# Write partial F row
writeLines("\\hline", file_conn)
writeLines(paste("Partial F:    &             &                       &", round(partial_F, 3),          "\\\\"), file_conn)
writeLines("\\hline", file_conn)
writeLines("\\hline", file_conn)
writeLines("\\\\", file_conn)

# End table and close file connection
writeLines("\\end{tabular}", file_conn)
close(file_conn)

  
rm("a", "b", "c_pipi", "c_deltadelta", "c_delpi", "crit", "cov_matrix", "pi", "delta", "n",  "k", "rss_full", "controls")